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Abstract 

We analyse the atomic state obtained by photo-dissociation of a molecular Bose- Einstein- 
condensate. This process is equivalent to down-conversion in quantum optics where it is 
responsible for squeezing of the field amplitudes. Monte Carlo simulations derived from the 
Positive P description of the system, and approximate equations for second moments of the 
atomic field operators are introduced to describe the early stages of the process - not amenable 
to the usual mean-field description. 

1 Introduction 

Like electromagnetic waves, atomic matter waves can to a large extent be manipulated in space. 
Interference phenomena have been observed, and a large number of high precision atom optics 
measurements have been carried out jl]] . One of the potential future applications of atomic Bose- 
Einstein condensates is as coherent matter wave sources for interferometric time and frequency 
standards, detection of inertial effects, and a host of related technological tasks. The spatial || 
and temporal || coherence of condensates has been verified, and coherent amplification of matter 
waves has been demonstrated Q] to establish the close analogy with laser and maser sources of 
light. 

For high precision purposes, however, it was realized long time ago, that so-called non-classical 
states of light may be more useful than classical field states H . It is therefore natural to consider 
the production of such 'non-classical' states of atoms as well. The term 'non-classical' is ambiguous 
here, since what is very non-classical for light, e.g. a number state, can seem perfectly classical 
for atoms, and we shall instead use the term quantum correlated states. Due to the collisional 
interaction between atoms in a Bose condensate, this system contains a non-linearity, in equivalence 
with the Kerr-effect for light, and this collisional interaction has already been observed Q in a 
matter wave analog of four-wave mixing in non-linear optics. There have been several proposals 
to utilize collisional effects to produce certain quantum correlated states, such as Schrodinger cat- 
like states 0, and to observe the effect of the non-linearity on the dynamics of the condensate 
||. It is clear from these studies that, even though the collisional interaction can be controlled 
experimentally ]j| it is not easy to control precisely the quantum features of interest; the problem 
is a genuine multi-mode one due to the spatial degrees of freedom, and one has to pay attention 
to the role of non-condensed atoms, as well. 

In this Paper, we suggest to follow instead the procedure of quantum optics for production 
of squeezed light: we suggest to implement an 'atomic OPO'. The optical parametric oscillator 
(OPO) is a device where light is down converted, so that a single pump photon at frequency 2uj 
is converted into two photons, each of frequency uj (degenerate case). The process is the inverse 
of second harmonic generation (SHG), and in practical experiments, one often sees a strong field 
of frequency lu, which is first frequency doubled in the SHG crystal, and the high frequency field 
is subsequently down-converted in a similar non-linear crystal, now working as an OPO. 
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The matter wave analogue of the SHG process is one in which the atoms are combined by photo- 
association into di-atomic molecules. This process has been analyzed theoretically |10[ |, and it was 
recently demonstrated experimentally, that part of an atomic condensate may be converted by 
stimulated Raman transitions into molecules this way (llj . We suggest to follow the analogy with 
the OPO quite directly, i.e., to apply laser fields to drive the Raman transition from the molecular 
state back to free atoms (photo-dissociation) , and we note, that the quantum correlations appear 
due to the fact that only even number states for the atomic component will be present in the 
sample, because atoms are created in pairs. 

The dynamics of coupled atomic and molecular degenerate gasses has been studied in some 
detail Q with particular focus on the spatial and temporal dependence of the mean field dynamics, 
e.g., the appearance of solitons We focus on the "quantum optics" aspects, i.e., the atom 

number distribution, and the consequences for atom counting experiments. 

The organization of this paper is as follows: In Sec. ||, we present our proposal, and we derive 
two theoretical methods, one approximate and one exact, to analyze the dynamics of the system. 
In Sec. ^ we present numerical results for various relevant quantities. In Sec. ^ we analyze the 
use of our quantum correlated atoms as squeezed input in two-state Rabi-oscillations with a large, 
"classical" , condensate. In Sec. [| we discuss the results and we briefly indicate some alternative 
ideas for the production and application of quantum correlated condensates. 



2 The model 

2.1 Trapped cold atoms in second quantization 

It is convenient to describe dilute gasses of cold atoms in a trap in the second quantized formalism. 
In this description the Hamiltonian for the system at low temperatures is usually taken to be: 



= I dr 



[ft{r)h^{r) + g -^{r)^{r)^{r)W)) . (1) 



The atomic field operators ip(f) and ip'(r) annihilate and create an atom at position r and they 
satisfy the equal time commutation relations 

tjj(r),ft{f)\ = S 3 (r-/) (2) 

The single particle Hamiltonian hi is the one appropriate for a single atom of mass m in the 
external trapping potential V ex t 

hi = -^V 2 + V ext (f). (3) 
2m 

The g term in Eq.([j]) describes the collisional two body interaction of the atoms. The simple 
contact form is an approximation appropriate for cold, dilute gases |l4| . The value of the strength 
parameter is g — A-Kh~a s /m where a s is the s-wave scattering length. 

In most experiments the trapping potential can be approximated well by a three dimensional 
harmonic oscillator characterized by three frequencies lu x , oj v and lu z . If one of these frequencies 
is significantly smaller than the other two the atoms will form an elongated cloud and we describe 
such an effectively one dimensional system by a single frequency to |L5| . 

We shall work with dimcnsionlcss equations and choose to measure time in units of u) , lengths 
in units of ao = y/h/uim, and energy in units of ftw. We are then left with 

H Q = Jdx {ft (x)hitp(x) + ^ft{x)$(x)$(x)${x)} (4) 

where 

1 92 1 2 

and where g — 4na s /ao is a dimensionally correct, one dimensional interaction strength [fusf . 
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2.2 Photodissociation from a molecular condensate 

Starting from a molecular condensate (prepared by photo-association or by other means) we can 
imagine to photo-dissociate or "down convert" the molecules to pairs of free atoms. A suitable 
Hamiltonian to describe photo-dissociation is 

This Hamiltonian clearly creates and annihilates atoms in pairs. A full description should also 
include the molecules such that each atomic pair creation would be accompanied by the annihi- 
lation of a molecule. Here we will assume the molecular condensate to be large and the number 
of molecules removed to be small. In that case it is reasonable to describe the molecular conden- 
sate by a time- independent c-number field. The function b in Eq. (|6j) is the product of this field 
and of the coupling to a position dependent laser field which is also assumed to be classical. To 
represent the position dependence of the molecular condensate and of the laser field, the relative 
wavefunction of the pair of atoms, and the time dependence of the laser field, we use the ansatz 

b(x,x',t) = exp (-2<At exp ^ exp -- V2V 7 

2ira r a cm \ 2 of J \ 2 a z cm ) 

where a r <C cr cm = <xq. 



2.3 Operator equations of motion 

A successful tool for the description of a Bose condensate is the Gross-Pitaevskii equation (GPE). 
It can be obtained from the Heisenberg equation of motion for the atomic field ip(x,t) by taking 
average values and by replacing the mean value of an operator product by the product of mean 
values. In our case the Heisenberg equation reads: 

. dip{x t) = fld + l 2+ ^^ t ^ x ^~^ + f b(x,x?,t)ft(xr,t). (8) 
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The last term in this equation is due to the exchange of atom pairs with the molecular condensate 
via photo-dissociation. 

It is easy to see that using the average of Eq.(||) has shortcomings when we try to describe 
photo-dissociation: If we start from the atomic vacuum defined by 

0)|0) = (9) 

the incoupling term on the right hand side has vanishing average value. Therefore (?p(x,t)) will 
stay zero also at all later times and no useful information can be extracted. 

To gain knowledge of the state created we therefore proceed to study expressions quadratic in 
the field operators. To shorten notation we define 

R(x,y,t) = ft(x,t)i>(y,t) (10) 
S(x,y,t) = j>(x,t)j>(y,t). (11) 

For these operators we get the following Heisenberg equations of motion: 

. dR(x,y,t) (Id 2 1 a 1 d 2 1 2 \ . 



dt \2dx 2 2 2 dy 2 2' 

+9 (R(y, y, t) - R{x, x, tn R(x, y, t) 

dz h(z,y,t)S t {x,z,t)-b*(x,z,t)S{z,y,t)\ (12) 
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and 

1 dt - \ 2dx 2 + 2 X 2dy* + 2 V ) b[X > V > t] 

+g (s(x -y) + R{x, x, t) + R(y, y, t)j S{x, y, t) 

+ J dz ^b{z,y, t)R(z,x) + b(x, z, t)R(z, y) j- + b(x, y, t). (13) 

Note that b(x, y, t) now appears as an inhomogeneous source term in the 5* equation. This guaran- 
tees a nontrivial behaviour when we take averages. Note also that when taking averages we have 
a problem with the interaction terms which will couple the second order expectations to fourth 
order expectations. These fourth order terms have to be factorized in some approximate way to 
obtain a closed set of equations. 

2.4 c-number equations 
2.4.1 Exact equations for g = 



When .9 = 0, Eq.(12) and Eq.([L3|) can be reduced to two coupled linear equations for the moments 
R(x,y,t) = {R(x,y,t)) , S(x,y,t) = {S(x,y,t)). (14) 



They read: 



: dR{x,y,t) (Id 2 1 2 1 d 2 1 2 



1 dt = {2d^-2 X -2W + 2 y ) RiX ' y ' t] 

+ dz {b(z,y,t)S*(x,z,t) — b*(x,z,t)S(z,y,t)} (15) 



dS(x,y,t) fid 2 1 2 Id , ., •, . 



dt \ 2dx 2 2 2 dy 2 2" 

+ J dz {b(z,y,t)R(z,x) + b(x,z,t)R(z,y)} 

+b{x,y,t). (16) 

Moreover, R and S uniquely determine all higher order expectation values. This can be seen in 
a number of ways. One is to note that the Wigner distribution is a multi-dimensional gaussian 
distribution fully characterized by its second order moments [0, [l^]. Another follows from the 
observation that when g = 0, the Heisenberg equation of motion (^|) implies that Tp(x, t) can at all 
times be expressed as a linear combination of the initial values 

$(x,t) = Jdy {f(x, y,t)i>(y, 0)+g(x, y,t)ft(y, 0)}. (17) 

Eq. (|i"7|) , its hermitian conjugate and the fact that our system starts in the vacuum state then 
suggest the following scheme for calculation of any operator product at arbritrary t: Use the 
commutation relation (||) to move all ^(x, 0) to the right of any fo(y, 0) (normal ordering). Of all 
the terms produced in this process only the ones consisting entirely of c-numbers are nonzero as 
the vacuum expectation of any normal ordered product of operators vanishes in the vacuum state. 
To evalute the c-number terms we formally need to calculate integrals of products of the / and 
g functions of Eq.(|l^). It is, however, not difficult to see that these integrals factorize and that 
the factors are exactly the ones involved in calculating expectation values of products of only two 
field operators. The end result is that the average of any operator product is replaced by a sum 
of all possible factorizations into two-operator expectations 

(ft(xi)i>Hx2)Tp(xS(x 4 )) = (^( X1 )^(X 2 ))WX 3 )^(X 4 )) 

+ {^(X 1 )^(X4))(^{X2)^{X 3 )). (18) 
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This is a simple version of Wick's theorem |l!| . 



2.4.2 Approximate equations for g ^ 

When g ^ we have to include the interaction term of Eqs.( ^2l^3|) in Eqs. ([l5|Jl^ ). Unfortunately, 
the decomposition Eq.(|l7|) is no longer exact, and there is no simple way to reduce the mean 
values of four-operator products to products of R and S. Rather than the simple replacement, 
e.g., R(y,y)R(x,y) — ► R(y,y)R(x,y), we choose to apply the Wick prescription as this is correct 
to lowest order. We then get 

g(R{y,y,t)R{x,y,t) - R{x,x,t)R(x,y,t)\ -> g(S*(x,y,t)S(y,y,t) - S(x,y,t)S*(x,x,t)) 

+2g {R{y, y, t) - R(x, x, t)) R{x, y, t) (19) 

(A(x, x, t) + R{y, y, t) + 5{x - y)) S(x, y,t)} -> 2g {R{y, y, t) + R(x, x, tj) S{x, y, t) 

+9 (R(x ) y,t)S(x,x,t) + R*{x,y,t)S(y,y,t)) 
+gS(x-y)S(x,y,t) (20) 

When these expressions are inserted into Eqs.([l5 16), we arrive at the equations we want to 



solve numerically. We use a split-step approach where the kinetic energy is treated by a Fourier 
method. The remaining terms are dealt with by a fourth order Runge-Kutta scheme. In this 
one-dimensional problem the equations are quite manageable. 



2.5 The positive P pseudo-probability distribution 

Pseudo-probablity distributions (PPD's) are well-established tools in quantum mechanics. The 
most well known of the distributions are the Glauber-Sudarshan P-function and the Wigner func- 
tion but especially in quantum optics a number of other distributions have also been useful. 
Common to all the PPD's is that they provide the expectation values of properly ordered operator 
products as weighted c-number averages. The positive P distribution (lj, [l5|, (P+) that we 
will be using here gives the expection values of normally ordered products by replacing ip by a 
c-number function tpi and "if)' by a c-number function -02, e.g.: 

Note that P + is the joint distribution of two spatial functions and it is therefore an immensely 
complicated functional in general. It satisfies a multi-dimensional Fokker-Planck equation which, 
however, opens the door to a Monte Carlo sampling as we can translate the Fokker-Planck equation 
for the distribution to Langevin equations for stochastic realizations of ipi (x, t) and ip2{x, t). These 
equations resemble the GPE but they are coupled and theycontain noise terms. A derivation of 
the equations without incoupling is given in [p^t (see also [pl| ) and the inclusion of incoupling is 
straightforward. In the notation of stochastic differential equations the equations can be written 



idipi(x) = (hiipi(x) + gip2(x)ijji (x)ipi(x)) dt (22) 

+ J b(x,x')^ 2 {x')dx' + dWi{x) (23) 

-id^ 2 {x) = (hi<ip 2 (x) + gipi(x)il> 2 (x)ip2(x)) dt (24) 

+ J b*(x,x')i/ji(x')dx' + dW 2 (x) (25) 

where h\ is still defined in Eq.fjlj) and the noise terms are gaussian and given by 

(dW h2 (x,t)) = (26) 

(dW 1 (x 1 t)dW 2 (x , ,t')) = (27) 

(dWi(x,t)dWi(x',t')) = idt{b{x,x',t)+gi; 1 (x,t)6(x-x'))S{t-t') (28) 

(dW 2 {x,t)dW 2 (x',t')) = -idt{b*{x,x',t) + gipi(x,t)S(x - x'))S(t- t'). (29) 
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By numerically simulating Eqs.( j25| , |25| ) we are able to calculate expectation values of arbitrary, 
normally ordered field operator products with the only approximation that the results are subject 
to sampling errors due to the use of finite ensembles. We describe in the appendix our procedure 
to synthesize the noise dW in our simulations. 

The crucial drawback of the method is the well know sudden divergence in some of the unphys- 
ical moments of the P+ distribution [pTj, |22] ]. Unphysical moments exist because the translation 
from operator products to products of ip\ and "02 never involves ip* and ip2- This leaves some 
room for P + to behave badly and unfortunately it exploits this freedom. In the wavefunction 
realizations, some wavefunctions diverge or they make very large excursions that are difficult to 
follow numerically and which make a devastating impact on the sampling error. The P+ Monte 
Carlo method is therefore limited to short times where only few atoms have been created and 
nonlinear effects are still small. This suits our purpose, since we are only interested in short time 
dynamics, and we shall trust the result produced by the Langevin equations as long as none of the 
wavefunctions in the ensemble have escaped the region where we have confidence in our integration 
algorithm. 



3 Results 

In this section we show results for some of the quantities of interest that we are able to calculate 
in our model. Although the main new feature lie in the quantum correlations we first show a 
very classical quantity, namely the density profile. We then proceed to look at the eigenvalues of 
the one-particle density matrix. The largest of these eigenvalues defines the condensate fraction 
and the corresponding eigenvector is the condensate wavefunction. Finally we turn to a two-body 
quantity, the second order correlation function. 

3.1 The density profile and the number of atoms 

The atomic density is given by the diagonal elements of the one-body density operator in the 
position representation, that is 



p(x,x) = {ip\x)i>{x)) =R(x,x) = ip2(x)ipi(x), (30) 

where the overbar in the last expression denotes the average over many realizations of the stochastic 
ipi(x,t) and ip2(x,t). In Fig. [I] is shown a typical plot of this profile at tot = 2.4. It has the 
characteristic gaussian shape of the harmonic oscillator ground state and as we shall see in Sec. |3.2| 
a large fraction of the atoms indeed occupy a common wavefunction close to this state. The R&S 
equations (|l5|,^6|) with the interaction terms ([l9],^0|) give results in excellent agreement with the 
P+ simulations. 

The total number of photo-dissociated atoms is obtained as the trace of the one-body density- 
operator or, according to Eq. (pfj|) , simply as 

N = J dx p{x,x). (31) 

In Fig. H this number is shown as a function of time for g = and for g — 0.01. The agreement 
between the R&S equations and the P + method is seen to be quite good. 



3.2 The condensate fraction and wavefunction 

By diagonalizing the one-body density matrix we obtain an orthonormal basis of single-particle 
states. The eigenvalues correspond to the populations of these states and in the case of a conden- 
sate one of these eigenvalues dominates, i.e., most particles occupy the same state. A dynamical 
picture of the condensation process is the Bose-enhancement of the scattering into the most oc- 
cupied state. Here we expect a similar effect to take place. At first several states of the system 
are occupied by the atoms created. As the number of atoms grows the stimulated character of the 
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creation becomes more important and a mode competition results in one mode being preferentially 
occupied. 

In Fig.|| we show the condensate fraction, i.e. the ratio of the largest eigenvalue of the one- 
body density matrix to the sum of the eigenvalues for different values of the interaction strenght 
g. It is seen that as expected the condensate fraction is in general an increasing function of time. 
The effects of interactions are rather small at these low atom numbers. Note that unlike studies 
of stationary condensates at T = 0, where interactions are responsible for the breakdown of a 
simple product state ansatz for the system and the existence of atoms outside the condensate, our 
incoupling by itself produces atoms both in the condensate and outside the condensate. In fact, 
our calculations show that the second-largest eigenvalue accounts for most of the atoms which are 
not in the condensate. 

As for the condensate wavefunction we see an interesting phenomenon: Although the density 
profile associated with the condensed part of the one-body density matrix is close to that of the trap 
ground state, the condensate wavefunction is in fact not stationary. The atoms have condensed into 
a state more resembling a squeezed statefj] and if the incoupling is stopped the wavefunction widths 
show an oscillating behavior. In Fig. [| we show (x 2 ) and (p 2 ) of the condensate wavefunction as 
a function of time. We see that at cut = 2.4 when the incoupling is stopped, the wavefunction is 
too wide in momentum space as compared to the ground state of the trap ({p 2 } > 1/2). 

One way to avoid this oscillation is to apply (5-kick cooling [^3| to the system. This procedure is 
efficient if there is a linear correlation between position and momentum. In the original suggestion 
the correlation between position and momentum is brought about by free expansion, but an 
examination of the condensate wavefunction shows that we have a similar correlation here. The 
idea is to apply a tight, harmonic trapping potential for a short time interval. If this interval 
is so short that any changes in position can be ignored, the effect is simply a momentum kick 
also varying linearly with position. Ideally, this kick brings all the particles to rest. In Fig. ^ we 
demonstrate that the procedure is effective in our problem. The results are shown for g = 0, but 
a similar reduction is achieved for non- vanishing g. 

3.3 The second order correlation function g( 2 >(x,y) 

More detailed information about the quantum state of the system is desired and available, and a 
natural quantity to consider is the second order correlation function 

It measures the probability to find two atoms at positions x and y normalized to the single-particle 
densities at the two positions. For a thermal state g^ = 2 while for a coherent state g^ = 1. 

As g^ involves the expectation value of a product of four field operators we are faced with 
similar factorization problem as when we derived the R&S equations. Again we will resort to the 
Wick prescription although it should be realized that this is only exact for states obtained with 
.9 = 0. 

We have already evaluated the expectation in the enumerator of Eq. (^) in terms of R and S in 
the gaussian case. This was done in Eq.(|l^) and we get for the second order correlation function 

(2)/ s , , \R(x,y)\ 2 + \S(x,y)\ 2 

9 (> (x,y) = l + ^7 t^t \ • (33) 

R{x,x)R(y,y) 

In contrast with the R&S equations the P+ method has no problems handling expectation 
values like the enumerator of Eq.([52|), and g^ can be determined exactly up to sampling errors. 
At relatively short times and low atom numbers we have therefore an excellent tool to obtain exact 
results even for g =/= 0. 

1 This position- momentum squeezing should not be confused with the the atom-field squeezing discussed later. 
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In Fig. U we show a plot of </ 2 **(0, 0) as a function of time for various values of g. The central 
value slightly above 3 indicates a strong bunching effect where two atoms are more likely to be 
found close together than in a coherent or a thermal state. This result can be compared with the 



analytical expression for a single mode squeezed state, generated by the Hamiltonian f3 ^a 



9 -+2 

2 -J- a T 



(2) (tftfaa) 1 ,„.s 

9 = TTTXrq = 3 + 7TTT7 (34) 

(ata) 2 (ata) v ; 

In the figure we plot both results of the R&S equations using Eq.(^) and the exact P + results. 
Good agreement is found between the two approaches until uot = 2.5 and hereafter the R&S 
equations fail to capture a decrease in the value of g^ 2 \ This decrease indicates a threshold effect 
that we will discus in the next section. 



3.4 Threshold effect 

In the semi-classical treatments of the laser and of the parametric oscillator, one identifies a thresh- 
old in the stationary balance between gain and loss; the fields shift from fluctuations around zero 
to fluctuations around finite intensities p4| . Above threshold these optical systems have smaller 
relative fluctuations of the intensity, and it is natural to expect a similar threshold behaviour in 
our model. There is a seemingly important difference between our model and the optical systems 
in the fact that we do not have an explicit dissipative mechanism. It has been known for a long 
time, and it has been demonstrated explicitly for a large number of physical systems, however, 
that quite generically, the interactions in many-body systems lead to ergodicity of eigenstates and 
a dynamical relaxation without coupling to an external bath. For a recent review, see |25j| . 

The R&S equations with their underlying assumption of a gaussian Wigner distribution, cen- 
tered around vanishing atomic field, are clearly unable to describe correctly the system around and 
above threshold, but the -P+ simulations are exact, and the discrepancy between the two meth- 
ods is thus most likely explained by a threshold effect. To investigate more closely the threshold 
hypothesis, we show in Fig. || scatter plots of ip2(0)4>i(0) at tot = 2.4, 3.0 and 3.6 for a situation 
with g = 0.02. From photodetection theory we can deduce the following expression for the atom 
number distribution 

P n (t) = — exp I — / ip2{x,t)ipi{x,t)dx I. (35) 

This expression, however, exhausts the statistical precision of the P + method, since it involves 
higher moments of the simulated amplitudes, which yield larger and larger fluctuations. Instead, 
we heuristically present histograms in the figures of how the real parts, Re(V>2(0)?/>i(0)) are dis- 
tributed, and these histograms are in good qualitative agreement with our picture of a bifurcation 
of the solution when we reach threshold in the process. It is seen how the distribution at ujt — 2.4 
is strongly peaked at zero with an exponential tail along the the real axis. At u>t = 3.0 this tail 
extends to larger values, and a shoulder around 120 atoms/ao starts to appear, and at tot = 3.6 a 
second maximum has developed. This explains the lowering of g^ seen in Fig. 0. 



4 Application of a squeezed condensate 

In this section we analyse a possible application of the state created in our model. We show how 
its peculiar statistical properties can be utilized to produce precise measurements. Our suggestion 
is in direct analogy with the use of a squeezed vacuum in quantum optics experiments with beam- 
splitters. When a beam in such an experiment is incident on a 50/50 beam-splitter and is split 
in two, the analysis of the noise properties (quantum fluctuations) of the two daughter beams 
depends crucially on realizing that the incoming beam is not only split at the beam-splitter, it is 
actually mixed with vacuum coming in from the back-side of the mirror. By replacing this vacuum 
by a squeezed state we may control the statistical properties of the daughter beams. 
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The matter wave analogue of the beam-splitter could in our case be a laser pulse which is 
able to coherently change the internal state of atoms. Such a pulse can drive each atom into a 
superposition of two internal states. 

Suppose now that we have two internal states of the atoms, state a and state b, and that we 
apply a 7r/2-pulse. We then have in the Heisenberg picture 

i>a = ^= tya + VV) (36) 

$a -> ^t, = -J= (j>a - V^fc) • (37) 

The total number operators of atoms in state a and state b after the pulse are thus given by 
N' a = I dxi>i{x)i>' a {x) 



~N a + l -N h + l -j dx [^ a {x)Mx) + ^l(x)Mx)} (38) 



Nl = I dx^l(x)^' b (x) 
1 1 a 1 



= -N a + -N b --j dx {^ a (x)Mx) + ^ h (x)Mx)} . (39) 
We will now concentrate on the difference in the number of atoms in the two states, 

K - N' b = f dx {ti(x)Mx) + $(x)Mx)} . (40) 



We will also let tjj a initially describe a large condensate while ipb describes our photodissociated 
state. That is, we assume that the photodissociation is to the internal state b, while at the time 
of the 7r/2-pulse these atoms are overlapped with a large normal condensate in the internal state 
a. The large condensate is assumed to be in a coherent state with a definite phase , e.g. 

Mx)\) = e ie cj> a {x)\) (41) 
m(x) = e-^OrXI (42) 

with 4> a real. The mean number of atoms is given by 

N a = (N a ) = J dxcf>l(x). (43) 

Using Eq.(|Fj|) we find 

'K -N' b )=Q (44) 



and 

((N' a -N^\ =N a + N b + J J dxdy Mx)My)^[R b {x,y,) + e- 2W S b {x,y,)] . (45) 

Ordinary vacuum in the 6-state is the special case with S b = R b — N b = and we note that the 
typical imbalance of populations is ^/Var(A^ — Nf) = VWa- Now, we use the nontrivial state 
created by photodissociation as squeezed vacuum. We imagine to have experimental control over 
and choose this phase optimally in order to reduce Var(iV^ — N b ). 

In Fig. [j] we plot this minimum value of Var(iV„ — N^) /N a as a function of the time of production 
of the squeezed condensate. It is clearly seen how the noise is rather quickly suppressed almost 
perfectly. In the particular case shown iV Q was taken to be 10 3 and <j> a was of gaussian shape. The 
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time-dependent number of atoms in the 6-condensate can be approximately read out of Fig. 
The width of <p a was chosen to be the same as the equilibrium width of the 6-state atoms. It is 
amazing, but of course already well-known in quantum optics, that this supression can take place 
even though the number of atoms initially in the 6-state is very small compared to the average 
number of atoms in the large condensate in the a-state 

At times beyond tot = 2 the noise suppression is lost in the exact P+ simulations. We recall 
that for g = 0, the R&S equations are exact, and it is thus natural to ascribe the discrepancy 
between the two methods to the interactions, which, by analogy with the Kerr-effect in optics, 
cause a deformation of the gaussian state. In Fig. || we show a scatter plot of the amplitudes 
of the projections, ipi a , of 3000 P+ realizations of tpi on <fi a at tot — 2.4 for g = 0.00 and for 
g = 0.02. Such plots should be interpreted with care, as in general the pair (tpi,^) is needed 
to calculate all normal-ordered expectations. The Kerr-effect, however, influences the S'-function, 
and to calculate (ip 2 ) we only need to average ip 2 . In our case the relevant contributions to this 
average can thus be depicted by a scatter plot of ip\ a alone. 

We see in Fig. ||a a gaussian state represented by points which form an ellipse-like structure. 
If the distibution had no prefered direction in phase-space (ip 2 ) = ip\ would vanish, but this is 
clearly not the case for our squeezed state. In Fig. |^b the the interaction modifies the phase 
accumulation of points with large values of \ipi\ and deforms the distribution to an S-like shape, 
and the mean value of ip 2 is reduced. If the state is centered around a finite field amplitude, the 
intensity dependent phase shift transforms a circular distribution into a bean-shaped one, and in 
that case the Kerr effect actually produces squeezing [p6[ . 



5 Discussion 



In this paper we have studied the matter wave analogue of an optical parametric oscillator. We have 
developed a description based on the correlations of the atomic field and we have supplemented 
this by the P + method. When interactions are taken into account both treatments were limited 
to relatively short times and low numbers of atoms. We showed how the state can be used as 
squeezed vacuum to improve the statistics of experiments with larger condensates. 

On the theoretical side, it is of course interesting to remedy the breakdown of the R&S method 
as we know that for large condensates an even simpler description, the Gross-Pitaevskii equation, 
is very successful. To bridge the gap between the two regimes which are both describable in simple 
terms would be interesting and could be useful also in other scenarios where special quantum states 
of a BEC are created. In such work it is very likely that time dependent Monte Carlo methods, 
based for example on the positive P distribution, P +1 will play a significant role. 

Let us finish this paper with a brief mentioning of other proposals for the preparation of 
quantum correlated atomic condensates. As mentioned in the introduction, the collisional inter- 
action has been proposed as an agent for the preparation of quantum correlated states, such as 
Schrodinger cat-like states (t]]. The states prepared by our proposal are substantially less 'non- 
classical', but provided the existence of a molecular condensate, we believe that the experimental 
requirements of our proposal are more easy to meet. 

It was suggested some time ago to prepare spin-squeezed states of non-degenerate atomic gasses 
and of trapped ions p7| , and concrete proposals were made, involving coupling of the atoms via 
their interaction with a quantized field mode (^8), or their center-of-mass motion [^9|. More re- 
cently, it was suggested, and experimentally proven possible, to induce transitions between atomic 
states in an optically thick sample by absorption of non-classical light, and thereby to transfer 
quantum correlations from the light field to the atoms Ingredients of these proposals be- 

come even more powerful with cavities around the samples |5l]. Very promising, recent ideas 
are based on quantum non-demolition (QND) measurements of, e.g., quadratures or populations 
in the atomic sample by refraction of light beams. Such detection suffices to establish quantum 
correlations in the sample, even though the detection is done with classical field states [B2|]. One 
may simply use the state following the QND measurement together with the known classical out- 
come of the measurement to perform subsequent high-precision experiments, or one may consider 
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feed-back loops, in which the system is driven towards a specific quantum correlated state. 

All of these proposals are also relevant for degenerate gasses. The QND methods are probably 
closest to real implementation, since phase contrast methods, already used for imaging of conden- 
sates in many experiments, only have to be carried out with properly chosen parameters to lead 
to quantum correlated atomic states. But all methods are interesting, and they may illuminate 
various aspects of the condensate dynamics, as illustrated, e.g. by the threshold phenomenon and 
the 'Kerr-effect' suppression of squeezing in the proposal in this paper. 



Synthesis of correlated noise 



In order to numerically simulate Eqs.(|23|,|25|) we discretize time and space, and we synthesize 
noise terms, dWi^fenjU), that obey discretized versions of Eqs.(^6||29|). It is well-known how 
independent (pseudo) random numbers from different distributions can be created. For example 
a Gaussian distribution can be created starting from uniformly distributed numbers via the Box- 
Miiller method, i.e. we know how to produce {dUi t 2(x n ,ti)} so that 

(dU h2 (x n ,U)) = (46) 
(dU a (x n ,U)dUp(x m ,tj)) = 5 (47) 



We see that the correlation functions, Eqs.(E8|29), contain two terms: one from the interaction and 



one from the incoupling. These can be treated separately if we split the noise in two independent 
contributions 

dW lt2 {x n ,ti) = dWl 2 (x n ,ti) + dWl 2 (x n ,ti) (48) 

Due to the contact form of the interaction, the corresponding noise term poses no difficulties; we 
simply choose 



dWl 2 (x n ,t t ) =^ ±i9 ^ i)dt dUl 2 (x n ,U), (49) 

where dUf 2 is chosen with the properties ( |46] , p7| ). The incoupling term is created by multiplication 
and convolution of uncorrelated noise d\J\ 2 of the form (|4^^) with suitable Gaussian functions: 

dWl 2 (x n ,ti) =AAexp(±# i A)exp (~^) Y. dx dU U x n> ,U)exp (- ^^) (50) 

It turns out that choosing 



2 o 2 2 2g c 2 m of / ±idtB 

4a^ m - cr^ y dxy/narCrcTnCrb 

is sufficient to fulfill Eq.(p6|) with b given by Eq. (]?]). o> and Ub are rather small, and in practice 
the sum in Eq.(pOj) only needs to involve a few terms. 
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Figure 1: The density profile at time uit = 2.4 for interaction strenght g = 0.01 and incoupling 
parameters B — 3.0h 2 /m, A — ui, a cm = cto and ay = 0.2ao- The solid curve shows the results of 
the coupled R&S equations, the crosses are obtained with the P + simulations. 
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Figure 2: The number of photo-dissociated atoms as a function of time. The solid and dashed 
curves show the results of the R&S equations for g = 0.00 and g = 0.01 while the symbols + and 
x indicate the corresponding results of the P + Langevin equations. Incoupling as in Fig.[j]. 
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Figure 3: The condensate fraction, i.e., the ratio of the largest eigenvalue of the one-body density 
matrix to the total number of atoms (trace of the one-body density matrix). Solid and dashed 
curves shown the results of the R&S equations for g — 0.01 and for g = 0.02. The corresponding 
results of the P+ simulations are indicated by the symbols + and x . The lower sets of data with 
the same symbols show the results when the incoupling is stopped at u>t = 2.4. 
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Figure 4: The freely evolving expectation values of (x 2 ) (solid curve) and (p 2 ) (dashed curve) are 
compared to the results obtained after application of a 5-kick to match the wavefunction with the 
ground state in the trap. If we apply a properly chosen kick just as we stop the incoupling at 
cut = 2.4 the oscillations can be almost completely removed. 
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Figure 5: The second order correlation function g^(x,y) evaluated at y — x = as a function 
of time. The solid and dashed lines show the results of the R&S equations for g = (exact) 
and g = 0.02. The symbols + and x show the results of the P + simulations with g = 0.00 and 
g = 0.02. The photo-dissociation stops at u)t = 2.4. For g ^ 0, here is a strong discrepancy 
between the results of the R&S equations and the P+ at times beyond u>t = 2.5. The effect is 
clearly dependent on the interaction. The symbol □ show P + results obtained when g — 0.02, and 
when the photodissociation is not interrupted. 
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Figure 6: Scatter plots of ^2(0)^1(0) from the P + realizations for u)t = 2.4 (a), u)t = 3.0 (b), 
and tot — 3.6 (c). The average of this quantity is ^ (O)-0(O)) , the central density in the trap and 
higher order, normally ordered, moments (ifj jin ip n ) are similar averages (^2(0)^1(0))". For clarity, 
we show in the bottom of the plots histograms of Re(i/>2(0)V'i(0))- It is seen how the character of 
the state changes with time as a second maximum in the distribution develops at a value different 
from zero. 
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Figure 7: The minimum value of Var(A^ — N' b )/N a as a function of the time of application of 
the it/ 2- pulse. All data are for a situation where incoupling is stopped at tot = 2.4. Calculations 
were done with three different interaction strengths (g = 0.005, 0.01 and 0.02) and the almost 
overlapping lines (solid, dashed and dotted) show the results of R&S equations while the symbols 
(+, x and 0) show the results of the P + simulations. 
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Figure 8: Scatter plots of V'la, the overlap of V>i with the condensate mode function (j> a /\fNa- 
Plot (a) is for g = 0.00 while plot (b) is for g = 0.02. The plots demonstrate the deforming effect 
of the interactions on the gaussian state. This deformation reduces the mean of ij)\ a which is the 
decisive factor for the efficiency of the state as squeezed vacuum. 
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